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Abstract 

The Saffman- Taylor viscous fingering instability occurs when a less 
viscous fluid displaces a more viscous one between narrowly spaced 
parallel plates in a Hele-Shaw cell. Experiments in radial flow ge- 
ometry form fan- like patterns, in which fingers of different lengths 
compete, spread and split. Our weakly nonlinear analysis of the insta- 
bility predicts these phenomena, which are beyond the scope of linear 
stability theory. Finger competition arises through enhanced growth 
of sub-harmonic perturbations, while spreading and splitting occur 
through the growth of harmonic modes. Nonlinear mode-coupling en- 
hances the growth of these specific perturbations with appropriate 
relative phases, as we demonstrate through a symmetry analysis of 
the mode coupling equations. We contrast mode coupling in radial 
flow with rectangular flow geometry. 

PACS numbers: 47.20.Ma, 47.20.Gv, 47.54.+r, 68.10.-m 

Keywords: Hydrodynamic instability; Interfacial dynamics; Pattern 
formation; Weakly nonlinear analysis. 
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1 Introduction 

The Saff man- Taylor problem Q addresses motion of two viscous fluids in the 
narrow space between two plates, known as a Hele-Shaw cell. When a fluid of 
low viscosity displaces a fluid of higher viscosity, the interface between them 
becomes unstable and starts to deform. This problem is formally equiva- 
lent to a special case of a solidification problem. Experiments and theory 
focus on two basic flow geometries (i) rectangular and (ii) radial ||. For 
both geometries, the initial developments of the interface instability tracks 
the predictions of linear stability theory |lj . As the unstable modes of pertur- 
bation grow, they become coupled in a weakly nonlinear stage of evolution. 
Finally, the system evolves to a complicated late stage, characterized by for- 
mation of fingering structures, in which nonlinear effects dominate. 

This paper addresses linear and the weakly nonlinear stages of the evo- 
lution. Here we concentrate our attention on the radial flow geometry. We 
develop the mode coupling theory for the Saffman- Taylor problem. Our ap- 
proach is a complementary study to the purely linear investigations |IJ and 
to recent developments in the analytical || and numerical || treatments of 
the fully nonlinear time-evolving flow. These analytical works |J describe 
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the early and late stage dynamics of the interface in the small surface ten- 
sion limit, using conformal mapping techniques. Efficient numerical simu- 
lations H describe complex flow patterns, for both large and small surface 
tension, employing a boundary integral formulation of the equations of mo- 
tion. On the other hand, our analytical weakly nonlinear study investigates 
the intermediate stage between the purely linear and the fully nonlinear ones, 
focusing on the onset of the nonlinear effects. Our approach applies to any 
value of surface tension, and gives physical insight into mechanisms of pat- 
tern solution. In the context of solidification, both numerical [^] and weakly 
nonlinear || |9| results are known. 

In a separate publication [|Hj] we develop the weakly nonlinear theory for 
rectangular flow geometry. Although mathematically similar to the rectangu- 
lar flow problem, the radial geometry imposes two important new elements: 
(i) the interface is always asymmetric because of the distinction between the 
area inside and outside of the boundary between the fluids; (ii) the radial flow 
presents multiple stages of instability, with no true steady state attainable. 

Experimental and theoretical investigations are plentiful for radial flow 
geometry Hele-Shaw cells These experiments [|], D], |I2, |13[ inject the less 



viscous fluid at the center of the cell, which was previously filled with the more 
viscous fluid. The injection is usually performed at a constant flow rate or 
constant pressure. Initially, the interface (bubble) is circular. As the bubble 
develops, the interface undulates and fingers grow. As the size of the bubble 
increases, the fingers spread and their tips become blunt, creating a fan-like 
pattern. Next, the fingers start to split at the tip. Larger fingers shield 
smaller ones and the pattern develops asymmetrically, causing competition 
among fingers of different lengths. Spreading, splitting and competition are 
the three basic growth mechanisms of the viscous fingering process [EH . 



The linear stability of a circular interface evolving in a radial geometry 



flow Hele-Shaw cell was studied by Bataille |13|], Wilson |TJ|], and thereafter 
by Paterson ||. More recent experimental and theoretical work study the 
linear stage of the radial Hele-Shaw flow in a variety of situations: Sader, 



Chan and Hughes [16] examined the linear stability analysis for the case of 



a non-Newtonian displaced fluid; Carrillo et al. [17] considered radial flow in 



a rotating Hele-Shaw cell; Cardoso and Woods [|18|] investigated the linear 
stability analysis of an annulus of fluid, bound by two other fluids of different 
viscosities. Because of the complexity of the evolution, even at the level of 



linear analysis, this paper reviews certain results of the linear theory for the 
sake of comparison with the nonlinear analysis. 

The linear stability analysis applies only to early stages of the flow. More- 
over, it does not explain well the essential mechanisms of spreading, splitting 
and competition cited above. To analyze these effects, we develop the weakly 
nonlinear mode coupling theory of the radial flow in a Hele-Shaw cell. We 
focus on the onset of the nonlinear effects and their influence on the evolu- 
tion of the unstable interface. We employ an analytical approach known as a 
mode coupling theory. The paper is organized as follows: section II derives a 
differential equation describing the early nonlinear evolution of the interface 
modes. In section III, we interpret results obtained in section II and investi- 
gate both the linear and weakly nonlinear evolution of the system. Our chief 
conclusions are presented in this section. We show that second order mode 
coupling explains the main features of patterning formation by coupling fun- 
damental modes to harmonics and sub-harmonics. Symmetries of the radial 
flow geometry dictate the form of the mode coupling equations, influencing 
the shapes of growing interfaces. Section IV presents our final remarks. 
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2 Derivation of the mode coupling differen- 
tial equation 

Consider two immiscible, incompressible, viscous fluids. In a Hele-Shaw cell 
fluids flow in a narrow gap of thickness b, between two parallel plates (see 
figure 1). We assume that b is smaller than any other length scale in the prob- 
lem, and therefore the system is considered to be effectively two-dimensional. 
Denote the viscosities of the inner and outer fluids, respectively as rji and 
T)2- The flows are assumed to be irrotational, except at the interface. Inject 
fluid 1 into fluid 2 at constant flow rate Q, equal to the area covered per unit 
time. Between the two fluids there exists a surface tension a. 

During the flow, the interface has a perturbed shape described as 1Z = 
R + C(0>*)) where 9 represents the polar angle, and R = R(t) denotes the 
time dependent unperturbed radius 



Ro is the unperturbed radius at t — 0. The analytic model we seek predicts 
the evolution of the interface perturbation amplitude £(9,t). The early non- 
linear evolution of the interface obeys a second-order mode coupling equation. 




(1) 
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The linear evolution is most readily described in terms of the complex 
Fourier modes 

Ut) = tt - r C(e, t) exp (-in6)d9, (2) 

Z7T JO 

where n=0,±l, ±2, ... denotes the discrete azimuthal wave number. We re- 
main in Fourier space for our weakly nonlinear analysis. The net perturbation 
is represented by 

+00 

C(M= E Ut)ex V (m9). (3) 

n=— 00 

Note that we include the n = mode in our Fourier expansion (|3]). This 
is done to keep the area A = nR 2 = -nR^ + Qt of the perturbed shape 
independent of the perturbation £. Thus we require 



1 00 

C0 = -^E[KnW| 2 +IC-n(t)| 2 ]. (4) 



The constant area constraint is intrinsically a nonlinear concern and is not 
required in linear stability analysis. Since we are interested in the early 
nonlinear behavior of the system, our first task is to derive a differential 
equation for £ n , correct to second order. 

We will occasionally refer to the "rectangular geometry limit" . By this we 
mean the limit R — > 00 and Q — ► 00 such that Q/(2ttR) = Voo and n/R = k 



remain constant, where v^, is the flow velocity at infinity and k denotes the 
wave number of the disturbance. In this limit the interface evolution reverts 
to the simpler evolution of the rectangular flow geometry. For example, the 
constant area constraint becomes simple (Co = 0) in this limit. 
The relevant hydrodynamic equation is Darcy's law |l|, |3|] 

- ^r** - VPi = o, (5) 

where Vi = Vi(r, 6) andpj = Pi(r, 6) are, respectively, the velocity and pressure 
in fluids 2 = 1 and 2. At the interface, the pressure difference between the 
two fluids is given by [|l[ 

(Pi - V-i) \n = o («|| + kj.) \n- (6) 
The curvature in the direction perpendicular to the plates (k±) is nearly 



constant |20|, ^T[]. Because its gradient is nearly zero, it does not significantly 
affect the motion in our problem. 

Due to the radial geometry, we have a somewhat complicated expression 



for the in-plane interface curvature |15|, 119 



where r denotes the distance from the radial flow source. The sign convention 
for the curvature Kn is such that a circular interface has positive curvature. 
Keeping second order terms in the perturbation amplitude (, we rewrite the 
in-plane curvature 



R R 2 \ dO 2 R 3 



C 2 + ^(^l +2C : 



1 fdCV , ^d 2 C 



(8) 



2 \ 09 J ^09 2 

For future reference, note that only one term (l/R 2 )(d 2 (/d8 2 ) out of six 
terms in equation (H) survives in the rectangular geometry limit. 

Taking advantage of the irrotational and incompressible flow conditions, 
we define the velocity potential 0; in each of the fluids, where Vi = — V0«. 
The velocity potential satisfies Laplace's equation V 2 0i = 0. Combining the 
velocity potential with equations (||) and (|7|) for the pressure difference and 
Darcy's law (H), we write the equation of motion 



A o - o = ~ aK b ( 9 ) 



(<t>l 


TZ + 4>2 






Tl-4>2 


A 


V 2 I 


\ 2 J 



where 



V2 + Vi 

is the viscosity contrast and 



A = ^1 (10) 



a= #° (11) 
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Equation (|9]) drops an arbitrary constant of integration related to the defini- 



tion of velocity potentials. For our weakly nonlinear analysis we are interested 
in second order contributions in the perturbation amplitudes. Therefore, 
all the quantities in equation (||) are evaluated at the perturbed interface 
1Z = R + C(9,t), and not at the initial unperturbed interface position R as 
is usually done in linearized surface deformation problems. The nonlinear 
terms arise because of this important distinction. 

Now define Fourier expansions for the velocity potentials fa. Far from the 
interface the velocity field should approach the unperturbed steady flow with 
a circular interface of radius R. Thus for r — > and r — > +00 the velocity 
potentials fa approach 0°, the velocity potentials for radial flow 



Ci are independent of r and 9. They do not contribute to the velocity fields Vi, 
so we will not discuss them further. The general velocity potentials obeying 
these requirements are 




(12) 




(13) 
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and 

P2n{t) [ ' ' 

The algebraic dependence on r transforms, in the rectangular geometry limit, 



+ (7) exp(infl). (14) 



into exponential dependence. In order to calculate the mode coupling differ- 
ential equation for the system, we substitute expansions fll^D , ( |13D and (|TjJ) 
into the equation of motion (Bp, keep second order terms in the perturbation 
amplitudes, and Fourier transform them. 

We need additional relations expressing the velocity potentials in terms of 
the perturbation amplitudes. To find these, consider the kinematic boundary 
condition which states that the normal components of each fluid's velocity 
at the interface equals the velocity of the interface itself 2^| . In polar coor- 



dinates (r, 0) the kinematic boundary condition is written as 
OTZ fldTZd^A (d<j>i\ 



dt \r 2 80 80 I „ \ dr I „ 



(15) 



Expand the boundary condition equation ( |T5D to second order in ( and then 
Fourier transform. The constraint that area A be independent of the pertur- 
bation ( (equation (f|)) enforces the equality of the n = Fourier mode on 
each side of equation (|15|). Nontrivial identities are obtained for 
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Solving for fa n (t) consistently to second order in ( yields 



Hn(t) I I Cn n i i Cn 

\n\ 2nK n 



+ \ sgn{nn') - j CwCn-n' + -^=^ sgn(nn')( n >(n 

\ \ n \J Z7lK «V0 



(16) 



and 



4>2n{t) = -r-Xn + p. r Cn 
|n| Z7rit|n| 

+ H ( Sgn{nn') + y^y ) Cn'Cn-n' + J2 S 9 n ( nn> )(n'(n-n'- 

(17) 

The overdot denotes total time derivative. The sign function sgn(nn') = 1 if 
(nn') > and sgn(nn') = —1 if (nn') < 0. We can use relations ( |T6"D and (|T7j) 
to replace the velocity potentials fa in Darcy's law @ with the perturbation 
C and its time derivative (. Keeping only quadratic terms in the perturbation 
amplitude, and equating Fourier modes n on each side of Darcy's law, leads 
to the differential equation for perturbation amplitudes ( n . For n^O, 



Cn = A(n)Cn + [ F ( n , n ')(ri(n-n> + G(n, n')Cn'Cn-n' , (18) 
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where 



A(n) 



^(A|n|-l)-^|n|(n»-l) 



(19) 



and 



|n| / rl 

r [n,n ) = — < 

V ' 7 R \2nR 2 



— — sgn[nn 



a 
R* 



n 



1 - —(3n' + n) 



G(n,ri) = -{A\n\[l - sgn{nn')\ - 1} . 
R 



(20) 
(21) 



Equation ( |Tq) is the mode coupling equation of the Saffman- Taylor prob- 



lem in a radial flow geometry Hele-Shaw cell. It gives us the time evolution 
of the perturbation amplitudes ( n accurate to second order. The first term 
on the right-hand side of equation (TXSh reproduces the linear stability anal- 



ysis P, O, ITSI], where A(n) denotes the linear growth rate. The remaining 



terms represent second-order mode coupling in a radial flow geometry. They 
arise as a direct consequence of our weakly nonlinear analysis, which consid- 



ers the presence of a full spectrum of modes. 



It is interesting to compare the radial geometry mode coupling equation 



with the rectangular geometry result 10 . Several terms simplify. In partic- 



ular, 



A(n) — > \k\ Avoo — ak 2 



(22) 
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F(n,n')-»0 (23) 
G(n, ri) -> A|ife| [1 - sflmCAA;')] • (24) 

When and fc' are parallel there is no second order mode coupling in the 
rectangular geometry limit. Note that, in contrast to the radial geometry 
expression (^9[) , the rectangular geometry linear growth rate ( ^2"D is not time 
dependent. The approach of each term to its limit is 0(1/ R). In the following 
section we investigate the radial geometry equation ([TJ) in more detail. 



3 Discussion 

3.1 Linear evolution 

Considering the complexity of the flow, we begin by analyzing certain aspects 
of the linear theory. This analysis is needed for our subsequent weakly nonlin- 
ear investigation (section 3.2). First consider the purely linear contribution, 
which appears as the first term on the right hand side of equation (|T8|). 
Since R varies with time, the linear growth rate A(n) is time dependent as 
well. This implies that the actual relaxation or growth of mode n is not 
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proportional to the factor exp[A(n)t], but rather 



Cn(t) =Cn(0)eX P 







X(n)dt' 



(25) 



If Jq \{n)dt' > the disturbance grows, indicating instability. Inspecting 
equation ([I9D for X(n) we notice opposing effects of the viscosity difference 
between the fluids (destabilizing) and of the surface tension at the inter- 
face (stabilizing). Two other relevant facts can be extracted from the linear 
growth rate: (i) the existence of a series of critical radii at which the interface 
becomes unstable for a given mode n (defined by setting A(n) = 0), 

„ . , 27ia \n\(n 2 — 1) 

«■<"> = (26) 

(ii) the presence of a fastest growing mode n*, given by the closest inte- 
ger to the maximum of equation (^) with respect to n (defined by setting 
d\(n)/dn = 0), 



n r . 



Note that n* varies with time. In view of equation (EH) n* is not simply 



related to the number of fingers present, even in the early stages of pattern 
formation. Furthermore, in the nonlinear regime the subsequent tip-splitting 
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process and mode competition result in a final number of fingers which can 
differ from the number present in the linear regime. 

Inspired by a model developed by Cardoso and Woods [18| (their "model 
B"), we describe the behavior of the system assuming the presence of a con- 
stant low level of noise during the whole evolution of the interface. The 
sources of noise may come, for instance, from inhomogeneities on the surface 
of the Hele-Shaw cell, irregularities in the gap thickness b, or even from ther- 



mal or pressure fluctuations [23|. The predictions of this model are in qualita- 



tive agreement with experimental observations within the linear regime 
We perform a detailed investigation of this model in the linear regime and 
extend its range of applicability to the weakly nonlinear stage of evolution. 

Suppose that we begin with an initially circular interface that is steadily 
expanding. During the interface expansion each mode n is perturbed with 
a constant (in time) random complex amplitude Cn(0). This noise ampli- 
tude contains an n dependent random phase but its magnitude |Cn(0)| is 
independent of n by assumption. As the interface continues to expand, it 
progressively reaches critical radii R c (n) for n — 2 , 3 , ... . Once a particular 
R c (n) is reached, the perturbation amplitude ( n starts to vary with time. 
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Within this model, the first order (linear) solution of equation (|I8" ) can be 
written as 



C n (t) 



Cn(0) if R < R c {n) 



Cn(0){(4^) A|nhl exp[(A|n|-l)(^-l)]} 



if R> R c (n). 

(28) 



To further analyze the behavior of the linear solution equation ( Pq) and 
see the overall effect of the above mentioned properties, we plot the time 
evolution of the interface. It is convenient to rewrite the net perturbation (H) 
in terms of sine and cosine modes 

oo 

C(0, t) = Co + E M*) cos ( n9 ) + b n{t) sm(n0)] , (29) 

n=l 

where a n = £ n + C-n an d b n = i (£ n — C_ n ) are real-valued. Throughout 
this work we use the experimental parameters given in Paterson's classical 
experiment [Q. Paterson observed the rapid growth of fingers, as air ~ 0) 
was blown at a relatively high injection rate, Q = 9.3 cm 2 /s, into glycerine 
(r/ 2 m 5.21 g/(cm s)) in a radial source flow Hele-Shaw cell. The thickness of 
the cell b = 0.15 cm and the surface tension a = 63 dyne/cm. We take into 
account modes n ranging from n = 2 up to 20. We evolve from initial radius 
Rq = 0.05 cm. The noise amplitude |Cn(0)| = -Ro/500. Figure 2, depicts the 
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evolution of the interface, for a random choice of phases, up to time t = 30 
s. 

Linear evolution proceeds through a cascade of modes, with increasing 
participation of higher modes n as time progresses. The number of fingers is 
typically given by the mode n which has grown to largest amplitude. Vari- 
ability of finger lengths depends primarily on modes of smaller n. Modes 
of larger n can cause the tips of the fingers to split. For the early evolu- 
tion of the interface, patterns like those shown in figure 2 resemble well the 



experimental patterns found in the literature [|], [TTJ, [Tj| . 

In contrast to the experiments, some important features are not present 
in figure 2, revealing failure of the purely linear approach. Spreading of the 
fingers, which is observed in the experiments, is not clearly shown in figure 2. 
Experimentally, the tips of fingers tend to split. To observe finger-splitting 
using the linear theory, we would have to go to unrealisticly long times. 
Finger competition, also, is not as pronounced in the linear theory as in 
experiment. These last few remarks indicate that, despite being an adequate 
approximation for early times, a purely linear theory does not describe all the 
basic physical mechanisms (spreading, splitting and competition) involved in 



the pattern formation. The linear regime is valid only for a limited initial 
period beyond which nonlinear effects take over. 

3.2 Weakly nonlinear evolution 

The nonlinear effects considered here come from the contributions of the sec- 
ond term on the right hand side of equation fllBp. Let us begin by briefly 
describing the various pieces appearing in this second term. The function 
F(n,n'), which multiplies the term Cn'Cn-n'; presents two parts (see equa- 
tion (|20D ) having distinct origins: the first one comes from the unperturbed 
radial flow (equation (|12|)) evaluated at the perturbed interface TZ, while the 
second one arises from the second order terms in ( present in the curvature 
(see equation (§)). Neither term survives in the rectangular geometry limit. 
On the other hand, the function G(n,n') (see equation (|21|)) multiplies the 
term Cn'Cn-n', which couples the perturbed flow ( with the interface shape 
perturbation (. Only the 0(1/ R) part of G(n,n') depends on the curvature 
of the unperturbed interface. 

To visualize the implications of mode coupling, we solve equation (|TS|) 
to second order accuracy. If we substitute the linear solution given in equa- 
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tion ( |28| ) into the second-order terms on the right hand side of equation (|T8|), 
we obtain the differential equation 

L = \(n)( n + W(n,t), (30) 

where 

W(n, *)= £ [>(", n)C l (^ n , + G(n, n')( l j(^ n ] (31) 

acts as a driving force in the linearized equation of motion (|30|) . Despite 
the complicated form of the functions A(n) and W(n, t), equation (|30"D is a 



standard first order linear differential equation [24| with the solution 



Cn(0) if R < R c (n) 

Ut) = { (32) 

Cn n (t) {l + It w [f^] if R > Rein). 
Here t c (n) is the time required for the unperturbed growth to reach radius 

R c {n) and can be easily calculated from equation ([!]). This solution describes 
the weakly nonlinear evolution, where the dominant modes just become cou- 
pled by nonlinear effects. 

We use the second order solution ( p2|) to investigate the nonlinear cou- 
pling among various modes n. Once again, we use the sine and cosine mode 
representation for the net perturbation ((8,t). In figure 3, we plot the in- 
terface for a certain time (t = 30 s), considering the same random choice 
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of initial phases as was employed in figure 2, and coupling all modes with 
2 < n < 20. The solid curve represents the interface shape obtained from 
the mode coupling solution ([32"]) , and the dashed line represents the interface 
in the linear approximation (taken from figure 2). The two interfaces are 
plotted together to facilitate comparison between the linear and nonlinear 
approaches. As we can see, the nonlinear evolution leads to wider fingers 
and their tips become more blunt. These fingers spread and some of them 
start to bifurcate, by splitting at the tip. The onset of finger competition is 
also observed, with some fingers slightly enhanced and others diminished. 

Even though the differences between the two evolutions are not dramatic, 
we see that the nonlinear approach developed here leads to a more realistic 
description of the radial fingering patterns. The effects of finger spreading, 
finger tip-splitting and finger competition are all present. The fact that the 
second order pattern is quantitatively close to the linear pattern suggests that 
the nonlinearities are fairly accurately described up to t = 30 s by second 
order mode coupling theory, because the contributions of third and higher 
order terms should be smaller yet. 

One interesting point to be discussed is the evolution of a particular mode 
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and its interaction with its sub-harmonics and harmonics (figure 4). Here 
we concentrate our attention on a small number of modes. This situation 
is obviously idealized compared with interaction of all modes illustrated in 
figure 3, but it gives us good indications about the role of mode coupling 
during interface evolution. 

For convenience we write the explicit equations of motion in the sine and 
cosine representation of the interface. The even and odd functions a n cos(n9) 
and b n sm(n9) form representations (respectively, the identity A±, and the 
reflection antisymmetric representation A 2 ) of the rotation and reflection 
symmetry group C nv f2"5|| . The signs, and relative magnitudes, of a n and b n 
determine the phase of mode n. In real space, they determine the location of 
the n fingers on the growing perimeter. In addition to these one-dimensional 
representations, C nv has a set of two dimensional irreducible representations 
{E m } formed by the cosine and sine functions of mode m — 1, 2, • • -, up to 
the largest integer less than n/2. Modes m' > n/2 transform identically to 
modes m = n — m' < n/2. When n is even, m = m' = n/2 is an integer, and 
the representation E n / 2 decomposes into a pair of reflection symmetric and 
antisymmetric representations Bi and B 2 . 
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The transformation properties of modes under rotations and reflections 
guarantee invariance of their equations of motion, but constrain the forms 
of those equations. Symmetry dictates which products of modes appear in 
each equation of motion. The numerical values of the coefficients multiplying 
these products, however, can only be obtained from a full derivation of the 
equation of motion such as we presented in section ^| Symmetry issues have 
been addressed in studies of shape selection for Saffman- Taylor fingers in 



rectangular channels [^, |27] and crystal growth processes |28|, |29|, |30 . 

First, we consider the influence of a fundamental mode on the growth 
of its harmonic. We take n as the fundamental and 2n as the harmonic. 
Without loss of generality we may take the phase of the fundamental so that 
a n > and b n = 0. We also replace a n with X(n)a n in terms already of 
second order. The equations of motion become 

a 2n = \{2n)a 2n + - [F{2n, n) + X(n)G(2n, n)} a\ (33) 



Jin 



X(2n)b 2n . (34) 



To see how symmetries constrain the form of these equations, we employ 
the group C 2n v The fundamental a n transforms as representation Bi, while 
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the harmonic a 2n transforms as the identity representation A\. The product 
of B\ with itself equals the identity representation, so a term of the form 
a n appears in equation (|33|). In contrast, b 2n transforms as the reflection 
antisymmetric representation A 2 , which cannot be formed from a product of 
the representation B\ with any other representation present. Consequently 
a n cannot influence the growth of b 2n at second order. 

Note that (1/2) [F(2n,ri) + X(n)G(2n, n)] < 0. At second order, the re- 
sult is a driving term of order a\ forcing growth of a 2n < 0. With this 
particular phase of the harmonic forced by the dynamics, the n fingers of 
the fundamental mode n tend, first, to spread and blunt, and later, to split. 
See figure 4 for an illustration. The solid line exhibits the broadened and 
split fingers predicted by the nonlinear mode coupling, while the dashed line 
shows the linear theory for comparison. Had growth of a 2n > been favored, 
fingers would have sharpened instead of blunting. Mode b 2n , whose growth is 
uninfluenced by a n , skews the fingers of mode n. In the presence of a 2n < 0, 
the role of h 2n is to favor one of the two split fingers over the other. The driv- 
ing term in equation of motion fl33|) spontaneously generates the harmonic 
mode, even in the absence of random noise. This effect is clearly shown in 
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the numerical solution of a solidification model 

In the rectangular geometry limit the driving term vanishes. Finger split- 
ting thus depends on the curvature of the unperturbed interface. Finger 
spreading and splitting is a unique prediction of the second order mode cou- 
pling in the radial geometry. 

Next, consider the influence of a fundamental mode n, assuming n is 
even, on the growth of its sub-harmonic mode n/2. Again, without loss of 
generality, we may choose the phase for the fundamental to set a n > and 
b n = 0. The equations of motion become 

a„/ 2 = |A(n/2) + - 



an^b n /2. (36) 

As in our analysis of the growth of the harmonic (equations ( |33"D and fl34|)) 
the multiplication of representations dictates the form of the equations of 
motion. We employ the symmetry group C nv . Mode n forms the identity 
representation A\. Note that a n / 2 and b n /2 transform as representations B\ 
and B2. Had we taken n odd, there would be no sub-harmonic n/2. However, 
the cosine and sine modes m transform as two-dimensional representations 
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2' 2 



\{n/2)G[l-l 



l n/2 



(35) 



5n/2 



A(n/2) 



n n\ w / 72 n 

/•'(--,- + A(n/2)G 



E m for m < n/2 and E n _ m i for m! > n/2. Consider a pair of modes m 
and m' = n — m. The product of the fundamental mode n with mode 
m! transforms as mode m. Consequently mode n and m' together would 
influence the growth of mode m. The special case under consideration has n 
even, and m = m! — n/2. 

Note that (1/2) [F(-n/2, n/2) + X(n/2)G(n/2, -n/2)) > 0. The funda- 
mental thus enhances the growth of its even sub-harmonic while inhibiting 
growth of its odd sub-harmonic. The result is to favor perturbations a n / 2 
which break the n-fold rotational symmetry of the fundamental by alter- 
nately increasing and decreasing the length of each of the n fingers. This 
effect describes finger competition. In contrast, mode 6 n / 2 skews the fingers 
while varying the depths of the valleys between fingers. Whereas the equa- 



tions of motion for the harmonic (equations fl33|) and (p4|)) contain a driving 
term, the equations of motion for the sub-harmonic (equations fl35|) and (|36"D ) 
show that rotational symmetry breaking sub-harmonics must be introduced 
through noise or initial conditions. They cannot grow spontaneously. The 
sign of a n /2 is not determined by the equation of motion. Positive and nega- 
tive values of a n /2 determine which fingers grow at the expense of the others. 
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Finger competition with a n /2 > is evident in figure 4. Reflection-symmetry 
breaking perturbations b n / 2 are suppressed consistent with known properties 
of selected steady state patterns |26|, |27| . 



Finally, consider the influence of a mode on itself. Consider a n > 
and b n = 0. There is no 0(a^) term in the equation of motion, consistent 
with the demands of invariance of the radial geometry equations of 
motion. However, our discussion above shows that a fundamental mode 
n spontaneously generates its own harmonic a 2n < of order a? n . This 
harmonic couples back into the equation of motion for a n (replace n/2 — > n 
in equation (pop), resulting in an 0(afy diminution over the linear growth 
of a n . In real space, a growing finger broadens as it pushes into the more 
viscous fluid. The broadened finger becomes increasingly resistant to further 
growth. The third order coupling of a mode with itself was studied in related 
solidification models ||, |[. 

In order to reinforce the conclusions of the previous paragraphs, compare 
the time evolution of the fundamental, harmonic and sub-harmonic pertur- 
bation amplitudes following the linear (|28|) and nonlinear (^) solutions. In 
figure 5 we present the time evolution of the perturbation amplitudes a n for 
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the modes n=4, 8 and 16 as given by the simple linear superposition of these 
modes (dashed curves) and by our weakly nonlinear mode coupling theory 
(solid curves). We observe that the weakly nonlinear coupling dictates the 
enhanced growth of both sub-harmonic (04) and first harmonic (ai 6 ) modes 
while diminishing the growth of the fundamental (a 8 ). The sign of the har- 
monic (die) is dictated as well, going strongly negative although its initial 
value was positive. Therefore, the nonlinear effects naturally enhance finger 
competition and tendency to finger tip-splitting. 



4 Concluding remarks 

In this paper, we studied the onset of nonlinear effects in the radial vis- 
cous fingering in a Hele-Shaw cell. While reasonably accurate for the initial 
stages of the flow, linear stability theories do not explain well the origin of 
the essential mechanisms of finger spreading, finger tip-splitting and finger 
competition. We derived a second order mode coupling differential equation 
for the perturbation amplitudes. By considering the interaction of a small 
number of modes, we revealed the role of sub-harmonic and harmonic per- 
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turbations at the onset of pattern formation. Symmetries of the radial flow 
geometry constrain the form of the mode coupling equations. Mode cou- 
pling provides a mechanism for finger spreading, splitting and competition 
creating patterns consistent with those observed in experiments. Reflection- 
symmetric fingers are favored. Our analytical mode coupling approach links 
the simple initial linear behavior with the strongly nonlinear advanced stages 
of evolution. 

We compared the radial flow geometry with the simpler case of flow in 
rectangular channels flIU| . The rectangular geometry represents a special 
case of the radial flow, in which effects of curvature of the unperturbed 
interface vanish. Finger splitting, in the radial flow geometry, is driven by 
curvature of the unperturbed interface. Tendency towards finger splitting 
may be controlled by adjusting this interfacial curvature. Presumably similar 
control is possible by performing the flow in curved two-dimensional spaces. 
Such effectively two-dimensional curved spaces may be found between curved 
surfaces embedded in three dimensions, separated from each other by a small 
distance b. Positive spatial curvature (e.g., flow between concentric spheres) 
should inhibit finger splitting while negative curvature (e.g., flow between 
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saddle shape surfaces) should enhance it. 

Finally, the weakly nonlinear approach we developed in this paper and 



in reference []10[] will be useful to investigate the onset of nonlinear effects 



in hydrodynamic stability problems involving ferrofluids p2| , plfl and lipid 



domains [32 . 
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Figure Captions 



Figure 1: Schematic configuration of the radial source flow. The viscosities 
of inner and outer fluids are rji and r] 2 , respectively. Fluid 1 is injected into 
the Hele-Shaw cell, previously filled with fluid 2, with constant injection rate 
Q. The dashed line represents the unperturbed circular interface R and the 
solid undulated curve depicts the perturbed interface 1Z = R + ((9,t), where 
6 is the polar angle. The surface tension between the two fluids is given by 
a and b denotes the thickness of the Hele-Shaw cell. 

Figure 2: Time evolution of the interface between the fluids, as given by 
the linear theory, including modes 2 < n < 20. The initial perturbation 
amplitudes |Cn(0)|=-Ro/500 and R = 0.05 cm. The experimental parameters 
are given in the text and £=10, 20 and 30 s. 

Figure 3: Snapshot of the interfacial patterns for t=30 s in the cases of 
linear (dashed curve) and nonlinear (solid curve) evolutions. All physical 
parameters and initial conditions are the same as those used in figure 2. 
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Figure 4: Snapshot of the interface between the two fluids (£=30 s), for the 
case of three interacting cosine modes (a fundamental n — 8, its sub-harmonic 
n/2 = 4 and its harmonic 2n = 16). Initial conditions are a n (0)=3.5 x 10~ 4 
cm and Rq = 0.3 cm. Other parameters are as in figures 2 and 3. The dashed 
(solid) curve obeys the linear (nonlinear) theory. The nonlinear contributions 
coming from the harmonic results in wider fingers, which tend to bifurcate. 
Finger competition is favored by the influence of the sub-harmonic mode. 

Figure 5: Time evolution of the cosine perturbation amplitudes a n (n=4, 
8 and 16) for linear (dashed curves) and nonlinear (solid curves) cases. The 
initial conditions and physical parameters are the same as those used in figure 
4. 
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